home *** CD-ROM | disk | FTP | other *** search
/ QRZ! Ham Radio 8 / QRZ Ham Radio Callsign Database - Volume 8.iso / mac / files / ant_nec / somnec.tz / somnec / somnec.for < prev   
Text File  |  1994-10-07  |  19KB  |  738 lines
  1. C     PROGRAM SOMNEC(INPUT,OUTPUT,TAPE21)
  2. C
  3. C     PROGRAM TO GENERATE NEC INTERPOLATION GRIDS FOR FIELDS DUE TO
  4. C     GROUND.  FIELD COMPONENTS ARE COMPUTED BY NUMERICAL EVALUATION
  5. C     OF MODIFIED SOMMERFELD INTEGRALS.
  6. C
  7.       PROGRAM SOMNEC
  8. C
  9.       COMPLEX CK1,CK1SQ,ERV,EZV,ERH,EPH,AR1,AR2,AR3,EPSCF,CKSM,CT1,CT2,C
  10.      1T3,CL1,CL2,CON
  11.       COMMON /EVLCOM/ CKSM,CT1,CT2,CT3,CK1,CK1SQ,CK2,CK2SQ,TKMAG,TSMAG,C
  12.      1K1R,ZPH,RHO,JH
  13.       COMMON /GGRID/ AR1(11,10,4),AR2(17,5,4),AR3(9,8,4),EPSCF,DXA(3),DY
  14.      1A(3),XSA(3),YSA(3),NXA(3),NYA(3)
  15.       DIMENSION LCOMP(4)
  16.       CHARACTER*32 OTFILE
  17.       DATA NXA/11,17,9/,NYA/10,5,8/,XSA/0.,.2,.2/,YSA/0.,0.,.3490658504/
  18.       DATA DXA/.02,.05,.1/,DYA/.1745329252,.0872654626,.1745329252/
  19.       DATA LCOMP/3HERV,3HEZV,3HERH,3HEPH/
  20. C
  21. C     READ GROUND PARAMETERS - EPR = RELATIVE DIELECTRIC CONSTANT
  22. C                              SIG = CONDUCTIVITY (MHOS/M)
  23. C                              FMHZ = FREQUENCY (MHZ)
  24. C                              IPT = 1 TO PRINT GRIDS.  =0 OTHERWISE.
  25. C     IF SIG .LT. 0. THEN COMPLEX DIELECTRIC CONSTANT = EPR + J*SIG
  26. C     AND FMHZ IS NOT USED
  27. C
  28. C     READ 15, EPR,SIG,FMHZ,IPT
  29.       PRINT100
  30. 100   FORMAT(' Program to Calculate Ground Interpolation Grid')
  31.       PRINT101
  32. 101   FORMAT(' For NEC2 Using Sommerfeld-Norton Method')
  33.       PRINT102
  34. 102   FORMAT(' ')
  35.       PRINT103
  36. 103   FORMAT(' Enter Relative Dielectric Constant:')
  37.       ACCEPT *, EPR
  38.       PRINT104
  39. 104   FORMAT(' Enter Conductivity (Mhos/Meter):')
  40.       ACCEPT *, SIG
  41.       PRINT105
  42. 105   FORMAT(' Enter Frequency (MHz):')
  43.       ACCEPT *, FMHZ
  44.       PRINT106
  45. 106   FORMAT(' Enter 1 to Print Grids, 0 to Suppress Printing:')
  46.       ACCEPT *, IPT
  47.       PRINT107
  48. 107   FORMAT(' Enter Data Output Filename:')
  49.       ACCEPT 24, OTFILE
  50.       PRINT *, ' Relative Dielectric Constant = ', EPR
  51.       PRINT *, ' Conductivity (Mhos/Meter) = ', SIG
  52.       PRINT *, ' Frequency, MHz = ', FMHZ
  53.       PRINT *, ' Printing Flag = ', IPT
  54.       PRINT *, ' Data Output File Name = ', OTFILE
  55.       IF (SIG.LT.0) GO TO 1
  56.       WLAM=299.8/FMHZ
  57.       EPSCF=CMPLX(EPR,-SIG*WLAM*59.96)
  58.       GO TO 2
  59. 1     EPSCF=CMPLX(EPR,SIG)
  60. 2     TST=SECNDS (0.0)
  61.       CK2=6.283185308
  62.       CK2SQ=CK2*CK2
  63. C
  64. C     SOMMERFELD INTEGRAL EVALUATION USES EXP(-JWT), NEC USES EXP(+JWT),
  65. C     HENCE NEED CONJG(EPSCF).  CONJUGATE OF FIELDS OCCURS IN SUBROUTINE
  66. C     EVALUA.
  67. C
  68.       CK1SQ=CK2SQ*CONJG(EPSCF)
  69.       CK1=CSQRT(CK1SQ)
  70.       CK1R=REAL(CK1)
  71.       TKMAG=100.*CABS(CK1)
  72.       TSMAG=100.*CK1*CONJG(CK1)
  73.       CKSM=CK2SQ/(CK1SQ+CK2SQ)
  74.       CT1=.5*(CK1SQ-CK2SQ)
  75.       ERV=CK1SQ*CK1SQ
  76.       EZV=CK2SQ*CK2SQ
  77.       CT2=.125*(ERV-EZV)
  78.       ERV=ERV*CK1SQ
  79.       EZV=EZV*CK2SQ
  80.       CT3=.0625*(ERV-EZV)
  81. C
  82. C     LOOP OVER 3 GRID REGIONS
  83. C
  84.       DO 6 K=1,3
  85.       NR=NXA(K)
  86.       NTH=NYA(K)
  87.       DR=DXA(K)
  88.       DTH=DYA(K)
  89.       R=XSA(K)-DR
  90.       IRS=1
  91.       IF (K.EQ.1) R=XSA(K)
  92.       IF (K.EQ.1) IRS=2
  93. C
  94. C     LOOP OVER R.  (R=SQRT(RHO**2 + (Z+H)**2))
  95. C
  96.       DO 6 IR=IRS,NR
  97.       R=R+DR
  98.       THET=YSA(K)-DTH
  99. C
  100. C     LOOP OVER THETA.  (THETA=ATAN((Z+H)/RHO))
  101. C
  102.       DO 6 ITH=1,NTH
  103.       THET=THET+DTH
  104.       RHO=R*COS(THET)
  105.       ZPH=R*SIN(THET)
  106.       IF (RHO.LT.1.E-7) RHO=1.E-8
  107.       IF (ZPH.LT.1.E-7) ZPH=0.
  108.       CALL EVLUA (ERV,EZV,ERH,EPH)
  109.       RK=CK2*R
  110.       CON=-(0.,4.77147)*R/CMPLX(COS(RK),-SIN(RK))
  111.       GO TO (3,4,5), K
  112. 3     AR1(IR,ITH,1)=ERV*CON
  113.       AR1(IR,ITH,2)=EZV*CON
  114.       AR1(IR,ITH,3)=ERH*CON
  115.       AR1(IR,ITH,4)=EPH*CON
  116.       GO TO 6
  117. 4     AR2(IR,ITH,1)=ERV*CON
  118.       AR2(IR,ITH,2)=EZV*CON
  119.       AR2(IR,ITH,3)=ERH*CON
  120.       AR2(IR,ITH,4)=EPH*CON
  121.       GO TO 6
  122. 5     AR3(IR,ITH,1)=ERV*CON
  123.       AR3(IR,ITH,2)=EZV*CON
  124.       AR3(IR,ITH,3)=ERH*CON
  125.       AR3(IR,ITH,4)=EPH*CON
  126. 6     CONTINUE
  127. C
  128. C     FILL GRID 1 FOR R EQUAL TO ZERO.
  129. C
  130.       CL2=-(0.,188.370)*(EPSCF-1.)/(EPSCF+1.)
  131.       CL1=CL2/(EPSCF+1.)
  132.       EZV=EPSCF*CL1
  133.       THET=-DTH
  134.       NTH=NYA(1)
  135.       DO 9 ITH=1,NTH
  136.       THET=THET+DTH
  137.       IF (ITH.EQ.NTH) GO TO 7
  138.       TFAC2=COS(THET)
  139.       TFAC1=(1.-SIN(THET))/TFAC2
  140.       TFAC2=TFAC1/TFAC2
  141.       ERV=EPSCF*CL1*TFAC1
  142.       ERH=CL1*(TFAC2-1.)+CL2
  143.       EPH=CL1*TFAC2-CL2
  144.       GO TO 8
  145. 7     ERV=0.
  146.       ERH=CL2-.5*CL1
  147.       EPH=-ERH
  148. 8     AR1(1,ITH,1)=ERV
  149.       AR1(1,ITH,2)=EZV
  150.       AR1(1,ITH,3)=ERH
  151. 9     AR1(1,ITH,4)=EPH
  152.       TIM=SECNDS (TST)
  153. C
  154. C     WRITE GRID ON TAPE21
  155. C
  156.       OPEN (UNIT=21,FILE=OTFILE,FORM='UNFORMATTED',STATUS='NEW',ERR=21)
  157.       WRITE (21) AR1,AR2,AR3,EPSCF,DXA,DYA,XSA,YSA,NXA,NYA
  158.       CLOSE (UNIT=21)
  159.       IF (IPT.EQ.0) GO TO 14
  160. C
  161. C     PRINT GRID
  162. C
  163.       PRINT 17, EPSCF
  164.       DO 13 K=1,3
  165.       NR=NXA(K)
  166.       NTH=NYA(K)
  167.       PRINT 18, K,XSA(K),DXA(K),NR,YSA(K),DYA(K),NTH
  168.       DO 13 L=1,4
  169.       PRINT 19, LCOMP(L)
  170.       DO 13 IR=1,NR
  171.       GO TO (10,11,12), K
  172. 10    PRINT 20, IR,(AR1(IR,ITH,L),ITH=1,NTH)
  173.       GO TO 13
  174. 11    PRINT 20, IR,(AR2(IR,ITH,L),ITH=1,NTH)
  175.       GO TO 13
  176. 12    PRINT 20, IR,(AR3(IR,ITH,L),ITH=1,NTH)
  177. 13    CONTINUE
  178. 14    CONTINUE
  179.       PRINT 16, TIM
  180.       GO TO 23
  181. 21    PRINT 22, OTFILE
  182. 23    STOP
  183. C
  184. 15    FORMAT (3E10.3,I5)
  185. 16    FORMAT (6H TIME=,E12.3,8H SECONDS)
  186. 17    FORMAT (30H NEC GROUND INTERPOLATION GRID,/,21H DIELECTRIC CONSTAN
  187.      1T=,2E12.5)
  188. 18    FORMAT (///,5H GRID,I2,/,4X,5HR(1)=,F7.4,4X,3HDR=,F7.4,4X,3HNR=,I3
  189.      1,/,9H THET(1)=,F7.4,3X,4HDTH=,F7.4,3X,4HNTH=,I3,//)
  190. 19    FORMAT (///1X,A3)
  191. 20    FORMAT (4H IR=,I3,/1X,(10(1PE12.5)))
  192. 22    FORMAT ('ERROR CREATING OUTPUT FILE = ',A)
  193. 24    FORMAT (A)
  194.       END
  195. C ***
  196. C
  197. C
  198.       SUBROUTINE BESSEL (Z,J0,J0P)
  199. C
  200. C     BESSEL EVALUATES THE ZERO-ORDER BESSEL FUNCTION AND ITS DERIVATIVE
  201. C     FOR COMPLEX ARGUMENT Z.
  202. C
  203.       COMPLEX J0,J0P,P0Z,P1Z,Q0Z,Q1Z,Z,ZI,ZI2,ZK,FJ,CZ,SZ,J0X,J0PX
  204.       DIMENSION M(101), A1(25), A2(25), FJX(2)
  205.       EQUIVALENCE (FJ,FJX)
  206.       DATA C3,P10,P20,Q10,Q20/.7978845608,.0703125,.1121520996,
  207.      1.125,.0732421875/
  208.       DATA P11,P21,Q11,Q21/.1171875,.1441955566,.375,.1025390625/
  209.       DATA POF,INIT/.7853981635,0/,FJX/0.,1./
  210.       IF (INIT.EQ.0) GO TO 5
  211. 1     ZMS=Z*CONJG(Z)
  212.       IF (ZMS.GT.1.E-12) GO TO 2
  213.       J0=(1.,0.)
  214.       J0P=-.5*Z
  215.       RETURN
  216. 2     IB=0
  217.       IF (ZMS.GT.37.21) GO TO 4
  218.       IF (ZMS.GT.36.) IB=1
  219. C     SERIES EXPANSION
  220.       IZ=1.+ZMS
  221.       MIZ=M(IZ)
  222.       J0=(1.,0.)
  223.       J0P=J0
  224.       ZK=J0
  225.       ZI=Z*Z
  226.       DO 3 K=1,MIZ
  227.       ZK=ZK*A1(K)*ZI
  228.       J0=J0+ZK
  229. 3     J0P=J0P+A2(K)*ZK
  230.       J0P=-.5*Z*J0P
  231.       IF (IB.EQ.0) RETURN
  232.       J0X=J0
  233.       J0PX=J0P
  234. C     ASYMPTOTIC EXPANSION
  235. 4     ZI=1./Z
  236.       ZI2=ZI*ZI
  237.       P0Z=1.+(P20*ZI2-P10)*ZI2
  238.       P1Z=1.+(P11-P21*ZI2)*ZI2
  239.       Q0Z=(Q20*ZI2-Q10)*ZI
  240.       Q1Z=(Q11-Q21*ZI2)*ZI
  241.       ZK=CEXP(FJ*(Z-POF))
  242.       ZI2=1./ZK
  243.       CZ=.5*(ZK+ZI2)
  244.       SZ=FJ*.5*(ZI2-ZK)
  245.       ZK=C3*CSQRT(ZI)
  246.       J0=ZK*(P0Z*CZ-Q0Z*SZ)
  247.       J0P=-ZK*(P1Z*SZ+Q1Z*CZ)
  248.       IF (IB.EQ.0) RETURN
  249.       ZMS=COS((SQRT(ZMS)-6.)*31.41592654)
  250.       J0=.5*(J0X*(1.+ZMS)+J0*(1.-ZMS))
  251.       J0P=.5*(J0PX*(1.+ZMS)+J0P*(1.-ZMS))
  252.       RETURN
  253. C     INITIALIZATION OF CONSTANTS
  254. 5     DO 6 K=1,25
  255.       A1(K)=-.25/(K*K)
  256. 6     A2(K)=1./(K+1.)
  257.       DO 8 I=1,101
  258.       TEST=1.
  259.       DO 7 K=1,24
  260.       INIT=K
  261.       TEST=-TEST*I*A1(K)
  262.       IF (TEST.LT.1.E-6) GO TO 8
  263. 7     CONTINUE
  264. 8     M(I)=INIT
  265.       GO TO 1
  266.       END
  267. C ***
  268. C
  269. C
  270.       SUBROUTINE EVLUA (ERV,EZV,ERH,EPH)
  271. C
  272. C     EVLUA CONTROLS THE INTEGRATION CONTOUR IN THE COMPLEX LAMBDA
  273. C     PLANE FOR EVALUATION OF THE SOMMERFELD INTEGRALS.
  274. C
  275.       COMPLEX ERV,EZV,ERH,EPH,A,B,CK1,CK1SQ,BK,SUM,DELTA,ANS,DELTA2,CP1,
  276.      1CP2,CP3,CKSM,CT1,CT2,CT3
  277.       COMMON /CNTOUR/ A,B
  278.       COMMON /EVLCOM/ CKSM,CT1,CT2,CT3,CK1,CK1SQ,CK2,CK2SQ,TKMAG,TSMAG,C
  279.      1K1R,ZPH,RHO,JH
  280.       DIMENSION SUM(6), ANS(6)
  281.       DATA PTP/.6283185308/
  282.       DEL=ZPH
  283.       IF (RHO.GT.DEL) DEL=RHO
  284.       IF (ZPH.LT.2.*RHO) GO TO 4
  285. C
  286. C     BESSEL FUNCTION FORM OF SOMMERFELD INTEGRALS
  287. C
  288.       JH=0
  289.       A=(0.,0.)
  290.       DEL=1./DEL
  291.       IF (DEL.LE.TKMAG) GO TO 2
  292.       B=CMPLX(.1*TKMAG,-.1*TKMAG)
  293.       CALL ROM1 (6,SUM,2)
  294.       A=B
  295.       B=CMPLX(DEL,-DEL)
  296.       CALL ROM1 (6,ANS,2)
  297.       DO 1 I=1,6
  298. 1     SUM(I)=SUM(I)+ANS(I)
  299.       GO TO 3
  300. 2     B=CMPLX(DEL,-DEL)
  301.       CALL ROM1 (6,SUM,2)
  302. 3     DELTA=PTP*DEL
  303.       CALL GSHANK (B,DELTA,ANS,6,SUM,0,B,B)
  304.       GO TO 10
  305. C
  306. C     HANKEL FUNCTION FORM OF SOMMERFELD INTEGRALS
  307. C
  308. 4     JH=1
  309.       CP1=CMPLX(0.,.4*CK2)
  310.       CP2=CMPLX(.6*CK2,-.2*CK2)
  311.       CP3=CMPLX(1.02*CK2,-.2*CK2)
  312.       A=CP1
  313.       B=CP2
  314.       CALL ROM1 (6,SUM,2)
  315.       A=CP2
  316.       B=CP3
  317.       CALL ROM1 (6,ANS,2)
  318.       DO 5 I=1,6
  319. 5     SUM(I)=-(SUM(I)+ANS(I))
  320. C     PATH FROM IMAGINARY AXIS TO -INFINITY
  321.       SLOPE=1000.
  322.       IF (ZPH.GT..001*RHO) SLOPE=RHO/ZPH
  323.       DEL=PTP/DEL
  324.       DELTA=CMPLX(-1.,SLOPE)*DEL/SQRT(1.+SLOPE*SLOPE)
  325.       DELTA2=-CONJG(DELTA)
  326.       CALL GSHANK (CP1,DELTA,ANS,6,SUM,0,BK,BK)
  327.       RMIS=RHO*(REAL(CK1)-CK2)
  328.       IF (RMIS.LT.2.*CK2) GO TO 8
  329.       IF (RHO.LT.1.E-10) GO TO 8
  330.       IF (ZPH.LT.1.E-10) GO TO 6
  331.       BK=CMPLX(-ZPH,RHO)*(CK1-CP3)
  332.       RMIS=-REAL(BK)/ABS(AIMAG(BK))
  333.       IF(RMIS.GT.4.*RHO/ZPH)GO TO 8
  334. C     INTEGRATE UP BETWEEN BRANCH CUTS, THEN TO + INFINITY
  335. 6     CP1=CK1-(.1,.2)
  336.       CP2=CP1+.2
  337.       BK=CMPLX(0.,DEL)
  338.       CALL GSHANK (CP1,BK,SUM,6,ANS,0,BK,BK)
  339.       A=CP1
  340.       B=CP2
  341.       CALL ROM1 (6,ANS,1)
  342.       DO 7 I=1,6
  343. 7     ANS(I)=ANS(I)-SUM(I)
  344.       CALL GSHANK (CP3,BK,SUM,6,ANS,0,BK,BK)
  345.       CALL GSHANK (CP2,DELTA2,ANS,6,SUM,0,BK,BK)
  346.       GO TO 10
  347. C     INTEGRATE BELOW BRANCH POINTS, THEN TO + INFINITY
  348. 8     DO 9 I=1,6
  349. 9     SUM(I)=-ANS(I)
  350.       RMIS=REAL(CK1)*1.01
  351.       IF (CK2+1..GT.RMIS) RMIS=CK2+1.
  352.       BK=CMPLX(RMIS,.99*AIMAG(CK1))
  353.       DELTA=BK-CP3
  354.       DELTA=DELTA*DEL/CABS(DELTA)
  355.       CALL GSHANK (CP3,DELTA,ANS,6,SUM,1,BK,DELTA2)
  356. 10    ANS(6)=ANS(6)*CK1
  357. C     CONJUGATE SINCE NEC USES EXP(+JWT)
  358.       ERV=CONJG(CK1SQ*ANS(3))
  359.       EZV=CONJG(CK1SQ*(ANS(2)+CK2SQ*ANS(5)))
  360.       ERH=CONJG(CK2SQ*(ANS(1)+ANS(6)))
  361.       EPH=-CONJG(CK2SQ*(ANS(4)+ANS(6)))
  362.       RETURN
  363.       END
  364. C ***
  365. C
  366. C
  367.       SUBROUTINE GSHANK (START,DELA,SUM,NANS,SEED,IBK,BK,DELB)
  368. C
  369. C     GSHANK INTEGRATES THE 6 SOMMERFELD INTEGRALS FROM START TO
  370. C     INFINITY (UNTIL CONVERGENCE) IN LAMBDA.  AT THE BREAK POINT, BK,
  371. C     THE STEP INCREMENT MAY BE CHANGED FROM DELA TO DELB.  SHANK'S
  372. C     ALGORITHM TO ACCELERATE CONVERGENCE OF A SLOWLY CONVERGING SERIES 
  373. C     IS USED
  374. C
  375.       COMPLEX START,DELA,SUM,SEED,BK,DELB,A,B,Q1,Q2,ANS1,ANS2,A1,A2,AS1,
  376.      1AS2,DEL,AA
  377.       COMMON /CNTOUR/ A,B
  378.       DIMENSION Q1(6,20), Q2(6,20), ANS1(6), ANS2(6), SUM(6), SEED(6)
  379.       DATA CRIT/1.E-4/,MAXH/20/
  380.       RBK=REAL(BK)
  381.       DEL=DELA
  382.       IBX=0
  383.       IF (IBK.EQ.0) IBX=1
  384.       DO 1 I=1,NANS
  385. 1     ANS2(I)=SEED(I)
  386.       B=START
  387. 2     DO 20 INT=1,MAXH
  388.       INX=INT
  389.       A=B
  390.       B=B+DEL
  391.       IF (IBX.EQ.0.AND.REAL(B).GE.RBK) GO TO 5
  392.       CALL ROM1 (NANS,SUM,2)
  393.       DO 3 I=1,NANS
  394. 3     ANS1(I)=ANS2(I)+SUM(I)
  395.       A=B
  396.       B=B+DEL
  397.       IF (IBX.EQ.0.AND.REAL(B).GE.RBK) GO TO 6
  398.       CALL ROM1 (NANS,SUM,2)
  399.       DO 4 I=1,NANS
  400. 4     ANS2(I)=ANS1(I)+SUM(I)
  401.       GO TO 11
  402. C     HIT BREAK POINT.  RESET SEED AND START OVER.
  403. 5     IBX=1
  404.       GO TO 7
  405. 6     IBX=2
  406. 7     B=BK
  407.       DEL=DELB
  408.       CALL ROM1 (NANS,SUM,2)
  409.       IF (IBX.EQ.2) GO TO 9
  410.       DO 8 I=1,NANS
  411. 8     ANS2(I)=ANS2(I)+SUM(I)
  412.       GO TO 2
  413. 9     DO 10 I=1,NANS
  414. 10    ANS2(I)=ANS1(I)+SUM(I)
  415.       GO TO 2
  416. 11    DEN=0.
  417.       DO 18 I=1,NANS
  418.       AS1=ANS1(I)
  419.       AS2=ANS2(I)
  420.       IF (INT.LT.2) GO TO 17
  421.       DO 16 J=2,INT
  422.       JM=J-1
  423.       AA=Q2(I,JM)
  424.       A1=Q1(I,JM)+AS1-2.*AA
  425.       IF (REAL(A1).EQ.0..AND.AIMAG(A1).EQ.0.) GO TO 12
  426.       A2=AA-Q1(I,JM)
  427.       A1=Q1(I,JM)-A2*A2/A1
  428.       GO TO 13
  429. 12    A1=Q1(I,JM)
  430. 13    A2=AA+AS2-2.*AS1
  431.       IF (REAL(A2).EQ.0..AND.AIMAG(A2).EQ.0.) GO TO 14
  432.       A2=AA-(AS1-AA)*(AS1-AA)/A2
  433.       GO TO 15
  434. 14    A2=AA
  435. 15    Q1(I,JM)=AS1
  436.       Q2(I,JM)=AS2
  437.       AS1=A1
  438. 16    AS2=A2
  439. 17    Q1(I,INT)=AS1
  440.       Q2(I,INT)=AS2
  441.       AMG=ABS(REAL(AS2))+ABS(AIMAG(AS2))
  442.       IF (AMG.GT.DEN) DEN=AMG
  443. 18    CONTINUE
  444.       DENM=1.E-3*DEN*CRIT
  445.       JM=INT-3
  446.       IF (JM.LT.1) JM=1
  447.       DO 19 J=JM,INT
  448.       DO 19 I=1,NANS
  449.       A1=Q2(I,J)
  450.       DEN=(ABS(REAL(A1))+ABS(AIMAG(A1)))*CRIT
  451.       IF (DEN.LT.DENM) DEN=DENM
  452.       A1=Q1(I,J)-A1
  453.       AMG=ABS(REAL(A1))+ABS(AIMAG(A1))
  454.       IF (AMG.GT.DEN) GO TO 20
  455. 19    CONTINUE
  456.       GO TO 22
  457. 20    CONTINUE
  458.       PRINT 24
  459.       DO 21 I=1,NANS
  460. 21    PRINT 25, Q1(I,INX),Q2(I,INX)
  461. 22    DO 23 I=1,NANS
  462. 23    SUM(I)=.5*(Q1(I,INX)+Q2(I,INX))
  463.       RETURN
  464. C
  465. 24    FORMAT (46H **** NO CONVERGENCE IN SUBROUTINE GSHANK ****)
  466. 25    FORMAT (10E12.5)
  467.       END
  468. C ***
  469. C
  470. C
  471.       SUBROUTINE HANKEL (Z,H0,H0P)
  472. C
  473. C     HANKEL EVALUATES HANKEL FUNCTION OF THE FIRST KIND, ORDER ZERO,
  474. C     AND ITS DERIVATIVE FOR COMPLEX ARGUMENT Z.
  475. C
  476.       COMPLEX CLOGZ,H0,H0P,J0,J0P,P0Z,P1Z,Q0Z,Q1Z,Y0,Y0P,Z,ZI,ZI2,ZK,FJ
  477.       DIMENSION M(101), A1(25), A2(25), A3(25), A4(25), FJX(2)
  478.       EQUIVALENCE (FJ,FJX)
  479.       DATA PI,GAMMA,C1,C2,C3,P10,P20/3.141592654,.5772156649,-.024578509
  480.      15,.3674669052,.7978845608,.0703125,.1121520996/
  481.       DATA Q10,Q20,P11,P21,Q11,Q21/.125,.0732421875,.1171875,.1441955566
  482.      1,.375,.1025390625/
  483.       DATA P0F,INIT/.7853981635,0/,FJX/0.,1./
  484.       IF (INIT.EQ.0) GO TO 5
  485. 1     ZMS=Z*CONJG(Z)
  486.       IF (ZMS.NE.0.) GO TO 2
  487.       PRINT 9
  488.       STOP
  489. 2     IB=0
  490.       IF (ZMS.GT.16.81) GO TO 4
  491.       IF (ZMS.GT.16.) IB=1
  492. C     SERIES EXPANSION
  493.       IZ=1.+ZMS
  494.       MIZ=M(IZ)
  495.       J0=(1.,0.)
  496.       J0P=J0
  497.       Y0=(0.,0.)
  498.       Y0P=Y0
  499.       ZK=J0
  500.       ZI=Z*Z
  501.       DO 3 K=1,MIZ
  502.       ZK=ZK*A1(K)*ZI
  503.       J0=J0+ZK
  504.       J0P=J0P+A2(K)*ZK
  505.       Y0=Y0+A3(K)*ZK
  506. 3     Y0P=Y0P+A4(K)*ZK
  507.       J0P=-.5*Z*J0P
  508.       CLOGZ=CLOG(.5*Z)
  509.       Y0=(2.*J0*CLOGZ-Y0)/PI+C2
  510.       Y0P=(2./Z+2.*J0P*CLOGZ+.5*Y0P*Z)/PI+C1*Z
  511.       H0=J0+FJ*Y0
  512.       H0P=J0P+FJ*Y0P
  513.       IF (IB.EQ.0) RETURN
  514.       Y0=H0
  515.       Y0P=H0P
  516. C     ASYMPTOTIC EXPANSION
  517. 4     ZI=1./Z
  518.       ZI2=ZI*ZI
  519.       P0Z=1.+(P20*ZI2-P10)*ZI2
  520.       P1Z=1.+(P11-P21*ZI2)*ZI2
  521.       Q0Z=(Q20*ZI2-Q10)*ZI
  522.       Q1Z=(Q11-Q21*ZI2)*ZI
  523.       ZK=CEXP(FJ*(Z-P0F))*CSQRT(ZI)*C3
  524.       H0=ZK*(P0Z+FJ*Q0Z)
  525.       H0P=FJ*ZK*(P1Z+FJ*Q1Z)
  526.       IF (IB.EQ.0) RETURN
  527.       ZMS=COS((SQRT(ZMS)-4.)*31.41592654)
  528.       H0=.5*(Y0*(1.+ZMS)+H0*(1.-ZMS))
  529.       H0P=.5*(Y0P*(1.+ZMS)+H0P*(1.-ZMS))
  530.       RETURN
  531. C     INITIALIZATION OF CONSTANTS
  532. 5     PSI=-GAMMA
  533.       DO 6 K=1,25
  534.       A1(K)=-.25/(K*K)
  535.       A2(K)=1./(K+1.)
  536.       PSI=PSI+1./K
  537.       A3(K)=PSI+PSI
  538. 6     A4(K)=(PSI+PSI+1./(K+1.))/(K+1.)
  539.       DO 8 I=1,101
  540.       TEST=1.
  541.       DO 7 K=1,24
  542.       INIT=K
  543.       TEST=-TEST*I*A1(K)
  544.       IF (TEST*A3(K).LT.1.E-6) GO TO 8
  545. 7     CONTINUE
  546. 8     M(I)=INIT
  547.       GO TO 1
  548. C
  549. 9     FORMAT (34H ERROR - HANKEL NOT VALID FOR Z=0.)
  550.       END
  551. C ***
  552. C
  553. C
  554.       SUBROUTINE LAMBDA (T,XLAM,DXLAM)
  555. C
  556. C     COMPUTE INTEGRATION PARAMETER XLAM=LAMBDA FROM PARAMETER T.
  557. C
  558.       COMPLEX A,B,XLAM,DXLAM
  559.       COMMON /CNTOUR/ A,B
  560.       DXLAM=B-A
  561.       XLAM=A+DXLAM*T
  562.       RETURN
  563.       END
  564. C ***
  565. C
  566. C
  567.       SUBROUTINE ROM1 (N,SUM,NX)
  568. C
  569. C     ROM1 INTEGRATES THE 6 SOMMERFELD INTEGRALS FROM A TO B IN LAMBDA.
  570. C     THE METHOD OF VARIABLE INTERVAL WIDTH ROMBERG INTEGRATION IS USED.
  571. C
  572.       COMPLEX A,B,SUM,G1,G2,G3,G4,G5,T00,T01,T10,T02,T11,T20
  573.       COMMON /CNTOUR/ A,B
  574.       DIMENSION SUM(6), G1(6), G2(6), G3(6), G4(6), G5(6), T01(6), T10(6
  575.      1), T20(6)
  576.       DATA NM,NTS,RX/131072,4,1.E-4/
  577.       LSTEP=0
  578.       Z=0.
  579.       ZE=1.
  580.       S=1.
  581.       EP=S/(1.E4*NM)
  582.       ZEND=ZE-EP
  583.       DO 1 I=1,N
  584. 1     SUM(I)=(0.,0.)
  585.       NS=NX
  586.       NT=0
  587.       CALL SAOA (Z,G1)
  588. 2     DZ=S/NS
  589.       IF (Z+DZ.LE.ZE) GO TO 3
  590.       DZ=ZE-Z
  591.       IF (DZ.LE.EP) GO TO 17
  592. 3     DZOT=DZ*.5
  593.       CALL SAOA (Z+DZOT,G3)
  594.       CALL SAOA (Z+DZ,G5)
  595. 4     NOGO=0
  596.       DO 5 I=1,N
  597.       T00=(G1(I)+G5(I))*DZOT
  598.       T01(I)=(T00+DZ*G3(I))*.5
  599.       T10(I)=(4.*T01(I)-T00)/3.
  600. C     TEST CONVERGENCE OF 3 POINT ROMBERG RESULT
  601.       CALL TEST (REAL(T01(I)),REAL(T10(I)),TR,AIMAG(T01(I)),AIMAG(T10(I)
  602.      1),TI,0.)
  603.       IF (TR.GT.RX.OR.TI.GT.RX) NOGO=1
  604. 5     CONTINUE
  605.       IF (NOGO.NE.0) GO TO 7
  606.       DO 6 I=1,N
  607. 6     SUM(I)=SUM(I)+T10(I)
  608.       NT=NT+2
  609.       GO TO 11
  610. 7     CALL SAOA (Z+DZ*.25,G2)
  611.       CALL SAOA (Z+DZ*.75,G4)
  612.       NOGO=0
  613.       DO 8 I=1,N
  614.       T02=(T01(I)+DZOT*(G2(I)+G4(I)))*.5
  615.       T11=(4.*T02-T01(I))/3.
  616.       T20(I)=(16.*T11-T10(I))/15.
  617. C     TEST CONVERGENCE OF 5 POINT ROMBERG RESULT
  618.       CALL TEST (REAL(T11),REAL(T20(I)),TR,AIMAG(T11),AIMAG(T20(I)),TI,0
  619.      1.)
  620.       IF (TR.GT.RX.OR.TI.GT.RX) NOGO=1
  621. 8     CONTINUE
  622.       IF (NOGO.NE.0) GO TO 13
  623. 9     DO 10 I=1,N
  624. 10    SUM(I)=SUM(I)+T20(I)
  625.       NT=NT+1
  626. 11    Z=Z+DZ
  627.       IF (Z.GT.ZEND) GO TO 17
  628.       DO 12 I=1,N
  629. 12    G1(I)=G5(I)
  630.       IF (NT.LT.NTS.OR.NS.LE.NX) GO TO 2
  631.       NS=NS/2
  632.       NT=1
  633.       GO TO 2
  634. 13    NT=0
  635.       IF (NS.LT.NM) GO TO 15
  636.       IF (LSTEP.EQ.1) GO TO 9
  637.       LSTEP=1
  638.       CALL LAMBDA (Z,T00,T11)
  639.       PRINT 18, T00
  640.       PRINT 19, Z,DZ,A,B
  641.       DO 14 I=1,N
  642. 14    PRINT 19, G1(I),G2(I),G3(I),G4(I),G5(I)
  643.       GO TO 9
  644. 15    NS=NS*2
  645.       DZ=S/NS
  646.       DZOT=DZ*.5
  647.       DO 16 I=1,N
  648.       G5(I)=G3(I)
  649. 16    G3(I)=G2(I)
  650.       GO TO 4
  651. 17    CONTINUE
  652.       RETURN
  653. C
  654. 18    FORMAT (38H ROM1 -- STEP SIZE LIMITED AT LAMBDA =,2E12.5)
  655. 19    FORMAT (10E12.5)
  656.       END
  657. C ***
  658. C
  659. C
  660.       SUBROUTINE SAOA (T,ANS)
  661. C
  662. C     SAOA COMPUTES THE INTEGRAND FOR EACH OF THE 6
  663. C     SOMMERFELD INTEGRALS FOR SOURCE AND OBSERVER ABOVE GROUND
  664. C
  665.       COMPLEX ANS,XL,DXL,CGAM1,CGAM2,B0,B0P,COM,CK1,CK1SQ,CKSM,CT1,CT2,C
  666.      1T3,DGAM,DEN1,DEN2
  667.       COMMON /EVLCOM/ CKSM,CT1,CT2,CT3,CK1,CK1SQ,CK2,CK2SQ,TKMAG,TSMAG,C
  668.      1K1R,ZPH,RHO,JH
  669.       DIMENSION ANS(6)
  670.       CALL LAMBDA (T,XL,DXL)
  671.       IF (JH.GT.0) GO TO 1
  672. C     BESSEL FUNCTION FORM
  673.       CALL BESSEL (XL*RHO,B0,B0P)
  674.       B0=2.*B0
  675.       B0P=2.*B0P
  676.       CGAM1=CSQRT(XL*XL-CK1SQ)
  677.       CGAM2=CSQRT(XL*XL-CK2SQ)
  678.       IF (REAL(CGAM1).EQ.0.) CGAM1=CMPLX(0.,-ABS(AIMAG(CGAM1)))
  679.       IF (REAL(CGAM2).EQ.0.) CGAM2=CMPLX(0.,-ABS(AIMAG(CGAM2)))
  680.       GO TO 2
  681. C     HANKEL FUNCTION FORM
  682. 1     CALL HANKEL (XL*RHO,B0,B0P)
  683.       COM=XL-CK1
  684.       CGAM1=CSQRT(XL+CK1)*CSQRT(COM)
  685.       IF (REAL(COM).LT.0..AND.AIMAG(COM).GE.0.) CGAM1=-CGAM1
  686.       COM=XL-CK2
  687.       CGAM2=CSQRT(XL+CK2)*CSQRT(COM)
  688.       IF (REAL(COM).LT.0..AND.AIMAG(COM).GE.0.) CGAM2=-CGAM2
  689. 2     XLR=XL*CONJG(XL)
  690.       IF (XLR.LT.TSMAG) GO TO 3
  691.       IF (AIMAG(XL).LT.0.) GO TO 4
  692.       XLR=REAL(XL)
  693.       IF (XLR.LT.CK2) GO TO 5
  694.       IF (XLR.GT.CK1R) GO TO 4
  695. 3     DGAM=CGAM2-CGAM1
  696.       GO TO 7
  697. 4     SIGN=1.
  698.       GO TO 6
  699. 5     SIGN=-1.
  700. 6     DGAM=1./(XL*XL)
  701.       DGAM=SIGN*((CT3*DGAM+CT2)*DGAM+CT1)/XL
  702. 7     DEN2=CKSM*DGAM/(CGAM2*(CK1SQ*CGAM2+CK2SQ*CGAM1))
  703.       DEN1=1./(CGAM1+CGAM2)-CKSM/CGAM2
  704.       COM=DXL*XL*CEXP(-CGAM2*ZPH)
  705.       ANS(6)=COM*B0*DEN1/CK1
  706.       COM=COM*DEN2
  707.       IF (RHO.EQ.0.) GO TO 8
  708.       B0P=B0P/RHO
  709.       ANS(1)=-COM*XL*(B0P+B0*XL)
  710.       ANS(4)=COM*XL*B0P
  711.       GO TO 9
  712. 8     ANS(1)=-COM*XL*XL*.5
  713.       ANS(4)=ANS(1)
  714. 9     ANS(2)=COM*CGAM2*CGAM2*B0
  715.       ANS(3)=-ANS(4)*CGAM2*RHO
  716.       ANS(5)=COM*B0
  717.       RETURN
  718.       END
  719. C ***
  720. C
  721. C
  722.       SUBROUTINE TEST (F1R,F2R,TR,F1I,F2I,TI,DMIN)
  723. C
  724. C     TEST FOR CONVERGENCE IN NUMERICAL INTEGRATION
  725. C
  726.       DEN=ABS(F2R)
  727.       TR=ABS(F2I)
  728.       IF (DEN.LT.TR) DEN=TR
  729.       IF (DEN.LT.DMIN) DEN=DMIN
  730.       IF (DEN.LT.1.E-37) GO TO 1
  731.       TR=ABS((F1R-F2R)/DEN)
  732.       TI=ABS((F1I-F2I)/DEN)
  733.       RETURN
  734. 1     TR=0.
  735.       TI=0.
  736.       RETURN
  737.       END
  738.